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ABSTRACT 

Previous  discussions  of  a  computing  method  for 
solving  two-dimensional  hydrodynamic  problems  are 
here  amplified  and  extended.  Results  of  computations 
are  presented  for  problems  involving  shock  diffraction 
and  refraction,  surface  instabilities,  and  viscous  flow. 


-  3  - 


TABLE  OF  CONTENTS 


Page 

Introduction  7 

Chapter  I.  The  Calculation  Procedure  10 

A.  Problems  Involving  Cartesian  Coordinates  in  a 

Rigid  Rectangular  Box  10 

1.  Layout  and  Nomenclature  10 

2.  The  Pressures  11 

3.  Phase  I  of  a  Calculation  Cycle  14 

4.  Phase  n,  The  Transport  of  Material  17 

5.  Phase  m,  Functionals  of  Motion  20 

B.  Other  Boundary  Conditions  in  Cartesian  Coordinates  21 

1.  Periodic  Channel  21 

2.  Prescribed  Input  21 

3.  Continuative  Output  22 

4.  Moving  Mesh  22 

5.  Rigid  Obstacle  23 

6.  Applied  Pressure  23 

C.  Generalized  Problems  in  Cartesian  Coordinates  24 

D.  Two-Dimensional  Calculations  in  Cylindrical  Coordinates  26 

E.  Limitations  of  the  Method  29 

Chapter  n.  Shock-Wave  Refraction  at  a  Gaseous  Interface  31 

A.  Introduction  31 

B.  Regular  Refraction  with  Corner  Effects  34 

C.  Irregular  Refraction  36 


-  5  - 


Page 

Chapter  m.  Shock  Passage  Through  a  Discontinuously 

Enlarged  Channel  45 

A.  Introduction  45 

B.  The  Computations  46 

C.  Development  of  the  Flow  Patterns  in  Nitrogen  47 

D.  Functionals  of  Motion  50 

E.  Infinite  Shock  in  Helium  51 

Chapter  IV.  Interaction  of  a  Shock  with  a  Deformable  Object  61 

A.  Introduction  61 

B.  Configurations  of  the  Flow  Fields  62 

C.  Functionals  of  Motion  62 

Chapter  V.  Hypersonic  Shear  Flow  with  Perturbed  Interface  68 

A.  Introduction  68 

B.  The  Interaction  69 

Chapter  VI.  Taylor  Instability  79 

Chapter  vn.  Viscous-Flow  Calculations  87 

A.  Introduction  87 

B.  Couette  Flow  87 

C.  Poiseuille  Flow  89 

References  86 


-  6  - 


INTRODUCTION 


The  particle-in-cell  method  for  two-dimensional  hydrodynamic  calcula¬ 
tions  has  been  applied  with  various  degrees  of  success  to  a  rather  wide 

variety  of  problems  in  compressible-fluid  flow.  The  method  was  first  dis- 

1  2 

cussed  in  two  unpublished  reports  *  which  have  been  superseded  by  more 

3  4 

detailed  published  descriptions  of  the  method  and  its  characteristics.  ’  The 
most  complete  description  previously  could  be  found  in  Ref.  4.  That  dis¬ 
cussion  was  mainly  restricted  to  one-dimensional  procedures,  however,  and 
as  applied  to  two-dimensional  calculations  was  incomplete  and  should  now  be 
modified  somewhat,  ft  is,  therefore,  one  purpose  of  this  report  to  discuss 
the  presently  used  procedure  in  some  detail. 

The  accuracy  of  the  computing  method  has  been  tested  by  applying  it 

to  a  variety  of  problems  for  which  theoretical  or  experimental  solutions  were 

5  6  7 

available.  The  results  of  some  of  these  calculations  have  been  reported;  ’  ’ 
others  are  available  only  in  classified  literature.  Still  others,  more  recently 
obtained,  have  revealed  new  restrictions,  or  regions  of  applicability,  or  have 
produced  results  not  previously  derived  by  theoretical  methods.  It  is,  there¬ 
fore,  the  second  purpose  of  this  report  to  summarize  these  new  results. 

For  brevity  in  writing,  the  particle-in-cell  method  for  hydrodynamic 
calculations  has  been  abbreviated  the  PIC  method.  In  the  discussions  to 
follow,  it  will  be  assumed  that  the  reader  has  access  to  Ref.  4  so  that  most 
of  the  discussions  presented  there  will  not  be  repeated.  On  the  other  hand, 
the  outline  of  computing  procedure  in  Chapter  I  of  this  report  is  sufficiently 


-  7  - 


complete  in  itself  so  that  a  computing  code  could  be  based  on  it. 

Performance  of  a  calculation  by  the  PIC  method  resembles  the  per¬ 
formance  of  an  experiment.  In  preparation,  the  differential  equations  of 
motion  are  transformed  to  suitable  conservative  finite-difference  forms. 

These,  together  with  the  initial  and  boundary  conditions  for  a  specific  situa¬ 
tion,  are  given  to  the  electronic  computer  which,  in  turn,  develops  the  solu¬ 
tion  at  a  sequence  of  later  times  separated  by  small  time  increments.  There 
is  no  a  priori  assumption  of  a  model  for  the  flow  configuration;  the  develop¬ 
ment  of  shocks,  for  example,  occurs  automatically  where  required.  Thus, 
these  computations  are  quite  different  from  another  type  often  performed  by 
high-speed  computers,  in  which  a  complicated  set  of  equations  is  solved  very 
precisely,  often  for  the  purpose  of  normalizing  analytical  approximation  pro¬ 
cedures.  Precise  solutions,  however,  are  usually  possible  only  with  ordinary 
differential  equations.  In  contrast,  the  PIC  method  approach  for  solving  the 
partial  differential  equations  of  hydrodynamic  always  results  in  approximate 
solutions.  It  has  been  observed  but  not  proved  that  under  many  circumstances 
of  interest  the  approximations  are  good  and,  furthermore,  that  they  can  be 
improved  by  decreasing  the  sizes  of  the  finite-difference  zones. 

In  the  absence  of  analytical  justifications  of  the  PIC  methodology,  It 
has  been  necessary  to  examine  by  "trial  and  error"  its  applicability  to  nu¬ 
merous  problems  with  known  solutions.  Likewise,  it  has  been  necessary  to 
experiment  with  numerous  modifications  of  the  methodology  in  order  to  obtain 
maximum  accuracy  with  the  resolution  presently  obtainable  with  available 
computing  machines.  In  some  cases,  it  has  been  found  that  a  modification 
would  result  in  very  little  change  in  the  answer;  In  other  cases,  a  small 
change  could  produce  a  very  large  effect.  Some  of  these  results  are  dis¬ 
cussed  in  this  report. 

The  calculation  method  is  designed  for  use  with  a  large  high-speed 


computer.  All  calculations  described  in  this  report  were  performed  on  an 
IBM  Electronic  Data  Processing  Machine,  type  704,  with  32K  memory. 

The  results  discussed  in  each  chapter  were  assembled  mainly  by  the 
people  mentioned  on  the  title  page.  There  were,  however,  many  occasions 
when  techniques  discovered  in  the  preparation  of  one  computing  code  were 
applied  to  another,  so  that  contributions  from  all  of  the  authors  can  be  found 
throughout  the  report.  In  addition,  contributions  to  these  later  developments 
in  the  PIC  methodology  have  been  made  by  Martha  W.  Evans  and  Billy  D. 
Meixner.  Many  stimulating  discussions  with  them  have  resulted  in  significant 
contributions  to  this  report;  details  of  their  recent  studies  are  to  be  reported 
elsewhere. 


CHAPTER  I 


» 


THE  CALCULATION  PROCEDURE 


A.  Problems  Involving  Cartesian  Coordinates  In  a  Rigid  Rectangular  Box 

1.  Layout  and  Nomenclature.  Two  materials  are  confined  to  move  in 
a  two-dimensional  rectangular  box  whose  walls  are  rigid  and  allow  perfect 
slippage.  The  materials  are  nonviscous  and  nonconducting  of  heat;  each  has 
an  equation  of  state  which  relates  pressure,  p,  to  density,  p,  and  specific 
internal  energy,  I. 

The  box  is  oriented  with  one  comer  at  the  origin  and  with  the  edges 
along  the  x  and  y  axes.  It  is  subdivided  into  a  number  of  equal  rectangular 

cells  to  which  the  finite -difference 
equations  are  to  be  related.  The 
cells  have  dimensions  fix  and  6y, 
whose  ratio  is  not  necessarily  the 
same  as  the  ratio  of  lengths  of  the 
box  sides.  A  typical  layout  is  shown. 

Each  fluid  is  represented  by  a 
number  of  mass  points  called  "par¬ 
ticles,"  each  with  a  constant  mass; 
as  shown  in  the  figure,  these  are 
represented  by  dots  and  x's;  we  shall 
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call  the  materials  "dot  material"  and  "x  material,"  respectively.  In  this 

example,  all  dot  particles  have  the  same  mass,  m. ,  and  all  the  x  particles 

have  the  same  mass,  m  .  (For  calculations  in  cylindrical  coordinates  or 

x 

for  certain  situations  in  cartesian  coordinates,  it  is  more  convenient  to  have 
a  different  mass  for  every  particle.)  For  each  particle  there  are  stored  in 
the  machine  memory  its  x  and  y  coordinates.  These  are  changed  in  time, 
by  the  method  described  below,  to  give  a  representation  of  the  motion  of  the 
fluids  through  the  mesh  of  cells. 

Such  quantities  as  velocity,  density,  and  pressure  are  kept  in  the 
machine  memory  by  cell,  so  that,  for  example,  the  pressure  of  a  cell  is 
a  certain  average  of  the  pressure  throughout  the  cell.  (Further  discussion 
of  this  point  is  given  by  Bromberg  in  Appendix  n  of  Ref.  4,  where  there 
is  an  enlightening  alternative  derivation  of  the  PIC  method  equations.)  The 
cells  are  labeled  with  index  with  i  and  j  increasing  in  the  x  and  y 
directions,  respectively;  the  lower  left  cell  in  the  figure  is  cell  number 
The  pressure  for  cell  (])  is  pj,  while  the  average  pressure  along  the 
boundary  between  cells  Q)  and  is  p|+i;  analogous  symbols  are  used 

for  the  other  boundary  pressures.  The  density  in  a  cell  is  defined  to  be 
the  quotient  of  the  sum  of  the  masses  in  the  cell  divided  by  its  area. 

The  nomenclature  for  various  cellwise  quantities  is  shown  in  Table 

1-1. 

2.  The  Pressures.  The  two  equations  of  state  are  given  in  the  form 
P  =  f  (P,  I) 

P  =  fx<P.D 

where  f  and  f  are  appropriate  functions  for  the  dot  and  x  materials,  re¬ 
spectively.  For  a  cell  containing  only  dot  or  x  material,  these  equations 
become 
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TABLE  1-1 


IDENTITY  OF  NOMENCLATURE  FOR  CELL 


uj  =  x-direction  component  of  velocity 
v|  =  y-direction  component  of  velocity 


M  J  =  mass  of  dot  material 
•i 


M  ^  =  mass  of  x  material  ( s  XT  ^ 
Xi 


(=  *  |m.) 

(s  N*S) 


(!) 


Xi 

A 


( )* 


M  *  +  M  J 
Xi 

specific  internal  energy  of  dot  material 
specific  internal  energy  of  x  material 
pressure 

the  Phase  I  change  in  total  internal  energy 

total  energy  of  dot  material 

total  energy  of  x  material 

total  x-direction  momentum 

total  y-direction  momentum 

result  of  Phase  I  calculation  for  (  ) 
result  of  Phase  n  calculation  for  (  ) 


M 


6x6y 


p  =  f 
i  x| 


M 


6x5y  * 


Various  procedures  are  possible  for  the  determination  of  total  pressure 
in  a  mixed  cell.  One  of  these  is  based  on  the  requirement  of  pressure  con¬ 
tinuity  across  a  material  interface.  Assuming  that  the  fraction  of  a  cell  oc¬ 
cupied  by  dot  material  is  <7,  one  writes  the  two  equations,  from  which  cr  is 
to  be  eliminated, 


M 


a<5x<5y 


=  f 


M 


~i 


I* 

(1  -  ar)6x6y  ’  x 


(1) 


If  the  pressure  is  strictly  proportional  to  the  density  for  both  materials, 
then  the  result  is  the  same  as  that  from  adding  partial  pressures: 


■ f. 


M 


6x6y  ’ 


+  f 
x 


6xSy  ’ 


(2) 


If  the  equations  of  state  are  complicated,  it  may  be  convenient,  as  well 
as  sufficiently  accurate,  to  still  calculate  the  pressure  in  a  mixed  cell  by 
adding  the  partial  pressures  in  this  manner.  In  some  cases,  however,  the 
result  of  this  is  far  from  reasonable  and  the  use  of  Eq.  (1)  is  indicated. 

An  approximate  solution  of  Eq.  (1)  has  been  found  useful  on  several  occasions. 
A  value,  c^,  is  assumed  for  a  and  the  pressure  is  taken  to  be 
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pj  - 1-1 
Fi  2 


M 


<T^6x6y 


+  f 


M 


x. 

1 


(l-o^xfiy’  xij 


The  choice  of 


M 


i 


j  j 

M  +  RM  J 
xi 


has  sometimes  yielded  reasonable  results,  where  R  is  the  ratio  of  the  initial 
density  of  the  dot  material  to  that  of  the  x  material.  Thus,  the  value  of 
is  based  on  the  assumption  that  the  compression  of  each  of  the  two  mate¬ 
rials  is  in  the  same  ratio  as  their  initial  compressions. 

Various  iterative  procedures  are  possible  for  solving  Eq.  (1);  none  of 
these  has  been  used. 


3.  Phase  I  of  a  Calculation  Cycle.  In  the  machine  memory  there 
are  stored  all  the  results  of  the  previous-cycle  calculations  or  else  the 
initial  conditions  for  the  problem.  These  are  to  be  advanced  in  time  accord 
ing  to  a  finite-difference  approximation  to  the  differential  equations 


dp 

at 


8u  au  ,  8u 

p - *•  pu  —  +  pv  —  + 

^at  H  ax  ay 


§e 

dx 


=  0 


av  ,  av  L  av 
P8t  +pUdk  +  pV^  + 
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t 


* 


The  first  of  these  equations,  that  of  mass  conservation,  is  automatically 
satisfied  by  the  particle  model.  The  other  three  equations  are  treated  as 
follows:  hi  Phase  I,  the  contributions  to  the  time  derivatives  which  arise 
from  the  terms  involving  pressure  are  calculated.  The  particles  are  not 
moved  at  this  step;  thus  the  transport  terms  are  dropped,  and  the  equations 
become,  in  finite  space -difference  form, 


From  these  result  the  tentative  new  cellwise  velocities 


~J  _  J  syfit  J  _  J  ‘ 
i  i  Pi+h  Pi-£ 


^.1  1  6x6t 


>*  -  pH 
pi  pi 


The  quantity  ij  in  the  third  equation  is  not  defined  for  a  mixed  cell,  but  we 
may  write 


ar! 


5Q 


i  at 


l  -  (-L-) 

t  \6x6y)  6t 


wherein  <5q|  is  the  change  in  total  internal  energy  of  the  cell.  Its  calcula¬ 
tion  involves  use  of  the  old  and  new  velocities 


-  15  - 


/  1  /  1\^ 

j- 

(“])  -W) 

-  6y6t 

-  0“>L 


—  6x6t 


j+i 

(pv):  2  -  (pv) 


H 

i 


It  is  seen  that  the  values  of  pressures  and  velocities  at  cell  boundaries 
are  required.  The  usual  procedure  in  obtaining  these  has  been  as  follows: 
Between  two  ordinary  cells,  the  boundary  pressure  and  velocity  are  taken  as 
the  averages  of  those  quantities  in  the  two  adjacent  cells.  At  the  rigid  box 
boundary,  the  normal  velocity  component  is  zero,  while  the  tangential  does 
/  not  need  specification.  The  pressure  at  the  box  boundary  has  been  taken  as 
a  quadratic  extrapolation  of  the  pressures  in  the  two  adjacent  cells.  The 
j  quadratic  extrapolation  is  used  because  it  allows  the  requirement  to  be  im- 
'  posed  that  the  normal  derivative  of  the  pressure  vanish  at  the  boundary.  It 
j  is  nearly  as  satisfactory  to  use  for  the  boundary  pressure  that  of  the  adjacent 
^  interior  cell. 

At  the  boundary  between  an  ordinary  cell  and  an  interior  empty  one, 
the  pressure  is  chosen  to  be  zero  so  that  the  velocity  there  need  not  be 
specified.  Such  a  choice  ensures  energy  conservation. 

It  has  been  found  to  make  little  difference  in  the  final  results  of  a 
calculation  whether  one  uses,  in  the  energy  equation,  the  product  of  average 
pressure  and  average  velocity,  or  the  average  of  the  pressure-velocity  pro¬ 
duct;  most  of  our  calculations  have  used  the  former. 

1  ^  1  ^  j 

From  the  result  for  6Q.,  the  values  of  I  '  and  I  are  to  be  determined. 

1  *  i  xA 

For  an  unmixed  cell  with,  say,  dot  material  only, 


+ 


<5Qj 


M 


3 

•i 
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One  way  of  distributing  internal  energy  in  a  mixed  cell  was  described 
in  Ref.  4  (page  21).  The  method  more  recently  used  is  based  on  the  as¬ 
sumption  that  both  components  are  compressed  or  expanded  adiabatieally 
through  the  same  pressure  change.  For  example,  if  both  components  are 
represented  by  polytropic  equations  of  state 

pk  “  <\  - 

(where  k  means  .  or  x),  then 


If  the  equations  of  state  are  more  complicated,  then  the  distribution  may  be 
based  on  other  assumptions,  such  as  (1)  both  components  receive  the  same 
change  in  total  internal  energy,  or  (2)  both  components  receive  the  same 
change  in  specific  internal  energy.  Application  of  the  first  of  these  two 
procedures  yielded  good  results  in  one  trial  calculation,  whereas  the  results 
of  using  the  second  were  not  as  good— insufficient  energy  was  transferred 
across  the  boundary  from  a  heavy  material  to  a  light  one. 

4.  Phase  n,  The  Transport  of  Material.  By  the  end  of  Phase  I, 
there  are  stored  in  memory  ten  quantities  for  every  cell.  Table  1-2  shows 
these,  together  with  the  quantities  which  replace  them  during  the  sequence 
of  Phase  n  calculations. 
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TABLE  1-2 


SEQUENCE  OF  STORAGE  CHANGES  FOR  CE LL ^  j  j  DURING  PHASE  H 

Step  1  > 

Step  2  > 

Step  3  > 

Step  4  > 


Step  1.  The  results  of  the  Phase  I  calculations  are  transformed  into 
total  cellwise  energies  and  momenta: 


X  =  (id  +  Mju 
Y  a  (.11  +  M  Jv 

Step  2.  The  particles  are  moved.  The  coordinates  of  each  change 
according  to 

x1  =  x  +  ufit 

y’  =  y  +  vfit 

hi  some  calculations  (see,  for  example,  Ref.  7),  the  values  of  u  and  v  were 
simply  the  values  of  u  and  v  of  the  cell  containing  the  particle,  no  matter 
where  in  the  cell  the  particle  originated  its  motion  for  the  cycle.  The  re¬ 
sults  can  almost  always  be  improved,  however,  by  using  the  time-consuming 
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process  called  'Velocity  weighting;"  in  fact,  tests  have  indicated  that  the  in¬ 
crease  in  accuracy  thereby  realized  could  not  be  achieved  by  that  increase 
in  mesh  fineness  which  would  consume  equal  extra  time. 

In  the  velocity  weighting  procedure,  a  rectangle  of  cell  size  is  imagined 
to  be  located  about  each  particle,  the  particle  being  at  the  center.  Such  a 
rectangle  then  overlaps  four  adjacent  cells,  and  the  effective  velocity  for 
moving  the  particle  is  taken  as  a  weighted  average  of  the  four  cellwise 
velocities,  the  weightings  being  proportional  to  the  overlap  areas.  If  the 
surrounding  rectangle  lies  partly  in  an  empty  cell,  then  that  cell  may  be 
assumed  to  have  the  same  velocity  as  does  the  cell  with  the  particle  (or 
any  other  convenient  and  reasonable  velocity)  for  the  purpose  of  determining 
the  effective  velocity  for  movement.  If  the  surrounding  rectangle  lies  out¬ 
side  the  walls  of  the  computation  region,  then  the  fictitious  outside  cells 
may  be  given  either  reflected  velocities  or  the  same  velocities  as  in  the 
adjacent  interior  cells.  In  the  former  case,  it  may  be  shown  that  for 
properly  small  values  of  St,  no  particle  will  be  lost  from  the  system.  The 
procedure  is  less  desirable,  however,  as  it  can  lead  to  the  "boundary  catas¬ 
trophe"  discussed  in  Ref.  4  (page  17).  In  the  latter  case,  it  is  necessary 
to  reflect  the  particle  back  in;  the  particle  then  carries  a  change  in  momen¬ 
tum  as  though  entering  from  a  cell  with  reflected  velocity,  and  the  boundary 
catastrophe  is  avoided. 

When  a  particle  is  thus  moved,  it  may  be  found  to  remain  in  the  same 
cell  from  which  it  started.  In  this  case,  there  is  no  modification  to  any  of 
the  cellwise  quantities.  Some  of  the  particles,  however,  will  end  up  in  new 
cells;  in  these  cases,  cellwise  changes  are  necessary.  From  the  cell  which 
was  left,  the  particle  mass,  momentum,  and  energy  are  subtracted  and  these 
are  added  to  the  new  cell.  Thus,  through  step  2,  the  cellwise  values  of 
mass,  momentum,  and  energy  cumulate  to  their  final  values  for  the  cycle. 


Step  3.  The  final  velocities  for  the  cycle  are  computed 
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Step  4.  The  final  specific  internal  energies  for  the  cycle  are  computed 
E«  . 
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5.  Phase  m,  Functionals  of  Motion.  The  final  arrangement  of  storage 
after  the  sequence  of  calculations  is  such  as  to  allow  immediate  re-entry 
into  Phase  I  of  the  next  cycle.  Ordinarily,  however,  before  proceeding  to 
the  next  cycle,  it  Is  useful  to  compute  various  functionals  of  the  motion  such 
as  total  kinetic  and  internal  energy  for  each  material,  components  of  total 
momentum,  positions  of  centers  of  mass,  entropy,  and  numerous  other  quan¬ 
tities.  In  some  cases,  it  is  possible  to  compare  changes  of  these  quantities 
with  the  changes  calculated  by  summing  boundary  fluxes.  Thus,  in  the  example 
at  hand,  the  total  energy  of  the  system  should  be  rigorously  conserved. 
(Actually,  machine  roundoff  will  introduce  some  change  in  total  energy,  but 
the  relative  change  per  cycle  should  be  bounded  by  a  number  which  is  pre¬ 
dictable  from  a  knowledge  of  the  number  of  significant  figures  retained  by 
the  calculation.)  Likewise,  changes  in  the  momentum  components  should  be 
exactly  predictable  in  terms  of  the  sum  of  the  boundary  forces.  Computed 
by  machine,  such  checks  serve  to  indicate  machine  or  coding  errors  and 
have  proved  extremely  valuable  on  many  occasions. 


B.  Other  Boundary  Conditions  in  Cartesian  Coordinates 


The  simple  boundary  conditions  discussed  in  the  previous  section  are 
applicable  only  to  a  restricted  class  of  problems.  Various  modifications 
are  listed  below;  most  of  them  have  been  used  in  actual  calculations. 

1.  Periodic  Channel.  The  rectangular  computation  region  can  be  con¬ 
sidered  to  be  one  section  of  an  infinite  channel  with  walls  parallel  to,  say, 
the  x  axis.  It  is  assumed  that  all  properties  of  the  entire  flow  field  are 
periodic  along  the  channel,  the  period  being  the  width  of  the  computation 
region.  The  change  in  computing  procedure  is  slight.  For  example,  along 
the  right-hand  boundary,  the  cells  are  treated  just  like  interior  cells  with 
their  right-hand  neighbors  being  the  cells  along  the  left-hand  boundary.  Par¬ 
ticles  leaving  the  system  across  the  boundary  re-enter  from  the  left,  while 
those  which  go  out  the  left-hand  boundary  are  inserted  from  the  right.  The 
condition  of  periodicity  is  applied  in  such  a  way  that  the  positions  of  the 
vertical  boundaries  are  immaterial  as  long  as  they  are  separated  by  one 
period.  Such  a  system  is  completely  conservative  of  particles,  energy,  and 
horizontal  momentum.  This  type  of  system  was  used  in  performing  the  cal¬ 
culations  discussed  in  Chapters  V  and  VII. 

2.  Prescribed  Input.  Along  one  or  several  of  the  boundaries,  a  pre¬ 
scribed  input  of  fluid  can  be  inserted.  This  could,  for  example,  be  used  to 
represent  the  flow  conditions  behind  a  shock  which  has  entered  across  one 
of  the  boundaries  at  the  beginning  of  the  problem.  Consider  the  example  of 
input  along  the  left  boundary.  The  left-hand  cells  of  the  computation  region 
can  be  treated  as  interior  cells  with  their  left-hand  neighbors  being  con¬ 
sidered  to  possess  the  prescribed  conditions  of  the  input  flow.  Particles 
are  periodically  ’'created"  for  insertion  across  the  left  boundary;  there  is 
thus  a  slight  additional  bookkeeping  difficulty  with  regard  to  the  storage  of 
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particle  coordinates.  This  type  of  input  system  has  been  used  in  numerous 
calculations,  including  those  discussed  in  Chapters  n,  HI,  and  IV,  and  in 
Refs.  5  and  7. 

3.  Continuative  Output.  Whenever  the  input  boundary  condition  (para¬ 
graph  2  above)  is  used,  a  provision  for  output  at  some  other  boundary  is 
usually  required.  If  the  flow  out  of  that  boundary  is  supersonic,  then  the 
exact  manner  by  which  it  is  treated  is  of  little  importance.  We  have  always 
used  a  '‘continuative”  boundary  treatment  for  such  an  output  line.  Accordingly, 
the  boundary  cells  are  treated  as  interior,  being  bounded  on  the  outside  by 
cells  with  identically  the  same  properties,  at  any  instant,  as  their  adjacent 
interior  neighbors.  The  machine  memory  locations  for  storage  of  the  co¬ 
ordinates  of  lost  particles  are  made  available  for  incoming  particles,  so 

that  the  total  required  machine  memory  is  bounded,  and  a  calculation  can 
be  continued  indefinitely. 

4.  Moving  Mesh.  In  paragraphs  1,  2,  and  3  above,  the  computation 
region  is  at  rest  and  the  fluid  streams  by.  Alternately,  if  it  is  desirable  to 
study  some  feature  of  the  flow  moving  with  fluid  speed  (or  some  other  speed), 
it  could  be  possible  to  have  a  traveling  region  of  computation.  Suppose,  for 
example,  one  wished  to  follow  the  motion  of  a  shock  wave  during  its  passage 
down  a  channel  in,  say,  the  x  direction.  This  could  be  accomplished  as 
follows.  A  zone  of  several  cells  would  always  be  present  ahead  of  the  shock; 
whenever  the  shock  had  advanced  a  cell  width,  a  new  column  of  cells  would 
be  created  to  their  right  with  conditions  representing  the  initial  state  ahead 
of  the  shock.  At  the  same  time,  a  column  of  cells  downstream  would  be 
destroyed.  Conditions  ahead  of  the  shock  could  be  constant  or  could  vary 
with  space.  Boundary  conditions  at  the  downstream  boundary  could  be  con¬ 
tinuative.  No  calculations  have  been  reported  which  used  a  moving  compu¬ 
tation  region  in  this  manner. 


5.  Rigid  Obstacle.  A  rigid  obstacle  can  be  placed  within  the  compu¬ 
tation  region.  This  is  most  easily  accomplished  if  the  boundaries  of  the 
obstacle  follow  cell  boundaries.  The  treatment  is  then  exactly  the  same  as 
at  the  rigid  walls  of  the  computation  region.  Such  a  calculation  was  reported 
in  Ref.  5.  If  the  obstacle  boundary  is  curved  or  oblique  relative  to  the  cell 
orientation,  then  the  procedure  is  somewhat  more  complicated.  Numerous 
partial  cells  are  created.  The  finite -difference  equations  for  such  cells  can 
be  derived  from  the  integral  form  of  the  equations  of  motion  by  a  procedure 
like  that  used  by  Bromberg  (Ref.  4,  Appendix  II)  for  deriving  the  equations 
under  ordinary  circumstances.  The  results  of  calculations  for  flow  past  a 
circular  object  were  reported  in  Ref.  7. 

6.  Applied  Pressure.  An  applied  pressure,  prescribed  as  a  function 
of  space  and  time,  can  be  exerted  on  the  irregular  surface  of  a  fluid.  We 
have  accomplished  this  by  varying  the  empty  cell  treatment,  using  the  empty 
cells  to  signal  the  presence  of  an  applied  pressure.  First,  the  motion  of 
particles  is  subjected  to  an  additional  constraint:  If  the  motion  of  any  par¬ 
ticle  results  in  emptying  a  cell,  then  the  particle  is  not  allowed  to  move 
during  that  cycle.  An  exception  is  allowed  if  the  emptying  cell  is  adjacent 
to  one  which  is  already  empty.  Initially,  the  fluid  interior  has  no  empty 
cells;  the  region  of  applied  pressure  is  an  empty-cell  region  surrounding 
the  fluid.  In  most  of  our  problems  the  resulting  motion  is  compressive,  so 
that  the  constraint  by  which  no  interior  cells  may  empty  is  not  serious. 

The  boundary  of  every  fluid  cell  adjacent  to  an  empty  cell  is  given  the  ap¬ 
propriate  applied  pressure,  and  it  is  assumed  that  the  velocity  at  that  bound¬ 
ary  is  that  of  the  fluid  cell.  The  pressure  within  a  cell  next  to  an  empty 
one  is  calculated  using  normal  density  in  the  equation  of  state,  if  the  cell 
would  otherwise  have  subnormal  density.  In  all  other  cases  these  edge  cells 
are  treated  as  ordinary  interior  cells.  This  procedure  has  been  used  in  the 


calculations  reported  in  Chapter  VI  and  in  other  calculations  not  reported. 

An  interpretation  of  the  applied  pressure  boundary  condition  is  as 
follows.  The  "empty"  cells  behave  as  though  they  were  filled  with  a  gas 
whose  density  is  very  small  compared  with  that  of  the  adjacent  material, 
but  whose  temperature  is  very  high  in  such  a  way  that  the  pressure  is 
finite  (the  prescribed  value).  As  a  result  the  sound  speed  is  very  high  so 
that  the  pressure  remains  homogeneous. 

C.  Generalized  Problems  in  Cartesian  Coordinates 

If  there  is  an  applied  body  force,  or  if  the  fluids  are  viscous  and  con¬ 
ducting  of  heat,  then  the  equations  and  boundary  conditions  become  somewhat 
more  complicated,  but  the  basic  method  is  not  altered.  As  an  example,  con¬ 
sider  the  problem  of  determining  the  nonsteady  motion  of  a  polytropic  gas 
flowing  through  a  periodic  two-dimensional  channel  bounded  by  walls  parallel 
to  the  x  axis.  The  basic  equations  (in  addition  to  that  of  mass  conservation, 
which  is  still  identically  satisfied)  are 


where  the  additional  symbols  are 
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first  (ordinary)  viscosity  coefficient 
second  viscosity  coefficient 

X/n  (assumed  constant;  an  idealized  monatomic  gas  has 

A-!) 


y  s  ratio  of  specific  heats  (assumed  constant) 

3  Prandtl  number  (assumed  constant) 

B  3  y/P^  [a  relatively  simple  theoretical  model  gives 

B  =  |(9y  -  5)] 

g  s  acceleration  of  body  force  (exerted  in  the  x  direction  only) 


Transformation  of  the  equations  to  finite -difference  form  proceeds  as 
before.  For  Phase  I  of  the  calculation 


(Equation  continued) 
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Again,  each  quantity  at  a  cell  boundary  is  taken  to  be  an  appropriate 
average  between  the  two  adjacent  cells.  Perhaps  the  only  feature  of  the 
finite-difference  form  which  is  not  otherwise  obvious  is  the  manner  in  which 
the  acceleration  term  is  written  in  the  energy  equation: 


rather  than,  say, 
uj  g6t 

Choice  of  the  form  shown  is  based  on  the  requirement  of  rigorous  energy 
conservation  by  the  difference  equations. 

More  details  for  a  specific  application  of  these  equations  are  given  in 
the  discussion  of  results  in  Chapter  VII. 

D.  Two-Dimensional  Calculations  in  Cylindrical  Coordinates 

There  are  many  problems  which  are  characterized  by  independence  of 
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the  flow  field  from  the  angle  about  some  fixed  axis.  Such  cylindrical  prob¬ 
lems  with  axial  symmetry  can  be  treated  by  the  PIC  method  with  almost  no 
modification  of  the  procedure  already  outlined.  Consider  the  problem  of 
determining  the  motion  of  a  single  nonviscous,  nonconducting  fluid  through 
an  infinite  periodic  cylindrical  pipe  with  rigid  walls.  In  some  plane  through, 
and  parallel  to,  the  axis,  one  period  of  the  flow  field  will  be  a  rectangle, 
bounded  at  the  bottom  by  the  axis,  at  the  top  by  the  pipe,  and  at  the  left  and 
right  by  the  ends  of  the  period.  This  region  is  divided  into  rectangular  cells 
(actually  toroids  of  revolution),  and  the  fluid  is  again  represented  by  particles 
(actually  circles  about  the  axis).  In  such  problems  it  has  been  found  con¬ 
venient  to  assign  different  masses  to  the  particles,  each  particle  being  given 
a  fixed  mass  whose  value  is  proportional  to  the  original  radius  of  the  par¬ 
ticle,  so  that  the  particle  density  is  initially  proportional  to  the  true  density. 

The  differential  equations 

du  _  _  8p 
Pdt  “  ‘dr 


dv  _  9p 


dt 


dz 


'it  [' +  i(u 


2  .  2) 


+  v 


1  9pur 
r  9r 


become,  for  the  Phase  I  calculations, 
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where  u  and  v  are  the  r-  and  z-direction  velocity  components,  respectively, 
and  i  and  j  count  cells  in  those  respective  directions.  At  the  axis,  u  =  0, 
while  a  quadratic  extrapolation  of  pressure,  as  discussed  before,  is  appro¬ 
priate.  These  difference  equations  are  perfectly  conservative.  Moreover, 
they  tend  to  preserve  spherical  symmetry.  This  fact  is  not  in  contradiction 
to  the  statements  in  Ref.  4  (p.  27)  wherein  the  difference  equations  were 
written  with  differences  between  cell  centers,  rather  than  between  cell  bound¬ 
aries.  The  present  form  of  the  equations  also  has  the  advantage  of  avoiding 
the  peculiar  boundary  condition  at  the  axis  required  in  Ref.  4  [Eq.  (52),  p.  28]. 

The  particle  movement  of  Phase  II  proceeds  by  a  velocity  weighting 
proportional  to  areas,  just  as  in  Cartesian  coordinates,  rather  than  to  vol¬ 
umes.  A  test  of  the  latter  procedure  produced  unsatisfactory  results,  es¬ 
pecially  near  the  axis.  Again,  if  a  particle  does  not  cross  a  cell  boundary, 
then  there  is  no  change  to  the  cellwise  quantities;  if  it  does  cross,  then 
mass,  "momentum,*'  and  energy  are  added  and  subtracted  as  before.  (Here 
the  radial-direction  "momentum"  is  defined  as  the  product  of  the  mass  and 
the  radial  velocity.) 

Numerous  calculations  using  the  PIC  method  in  cylindrical  coordinates 
show  that  it  is  nearly  as  useful  for  such  problems  as  for  those  in  Cartesian 
coordinates,  even  when  there  may  be  large  radial-direction  motions  of  the 
fluid.  Results  of  one  series  of  cylindrical  problems  are  discussed  in  Chap¬ 
ter  VI. 
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E.  Limitations  of  the  Method 


The  PIC  method  has  been  found  useful  in  solving  a  wide  variety  of 
problems  concerning  the  dynamics  of  compressible  fluids.  Very  little  ana¬ 
lytical  work  has  been  accomplished  in  the  direction  of  proving  the  validity 
of  results,  so  that  considerable  experimentation  has  been  necessary.  The 
range  of  applicability  is  discussed  in  this  report  and  elsewhere;  it  is  also 
appropriate  again  to  emphasize  the  limitations. 

The  principal  limitation  of  the  PIC  method  arises  from  the  requirement 
that  the  fluid  speed  relative  to  the  computational  mesh  must  not  be  small 
compared  with  the  sound  speed.  An  exception  is  allowed  in  a  uniform  region, 
where  the  fluid  speed  may  be  zero.  Thus  it  is  not  possible  to  apply  the 
method  to  problems  in  incompressible-fluid  flow.  There  are  two  related 
reasons  for  this  restriction.  First,  interactions  within  the  fluid  are  propa¬ 
gated  only  from  cell  to  cell,  whereas  in  an  incompressible  fluid,  the  changes 
in  configuration  at  a  point  depend  upon  conditions  at  that  instant  throughout 
the  fluid.  Second,  as  the  fluid  speed  decreases,  the  "effective  viscosity"  due 
to  the  dissipative  procedure  used  in  cell  crossings  (Ref.  4,  p.  16)  decreases 
to  zero.  The  difference  equations,  in  this  limit,  can  be  shown  to  be  uncon¬ 
ditionally  unstable.  Thus,  in  a  region  of  "perturbed  stagnation,"  fluctuations 
of  the  field  variables  grow  until  the  velocities  are  large  enough  to  produce 
sufficient  dissipation.  A  further  discussion  of  this  matter  is  given  in  Ref.  6, 
page  10.  This  second  difficulty  can  be  relieved  somewhat  by  the  incorpora¬ 
tion  of  artificial  dissipative  terms  into  the  difference  equations.  Usually 
such  terms  also  have  somewhat  undesirable  results,  such  as  the  increased 
smearing  of  discontinuities.  A  discussion  of  these  extra  terms  is  in  prepa- 

g 

ration  by  Longley. 

Another  significant  limitation  of  the  PIC  method  results  from  the  in¬ 
ability  of  the  fixed  mesh  of  cells  to  resolve  features  of  the  flow  field  which 
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are  small  in  size  compared  to  the  over-all  region  to  be  studied.  In  some 
cases  the  limitation  can  be  overcome  by  the  creation  and  destruction  of 
cells,  so  that  computational  mesh  is  present  only  where  needed.  It  would 
also  be  possible  to  have  fine  zones  at  some  localities  and  coarse  zones  at 
others,  but  such  a  procedure  has  not  yet  been  tried. 
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CHAPTER  H 

SHOCK-WAVE  REFRACTION  AT  A  GASEOUS  INTERFACE 

A.  Introduction 

9 

In  a  recent  paper,  Jahn  presented  the  results  of  a  set  of  experiments 
designed  to  study  the  regular  and  irregular  refraction  patterns  arising  from 
the  interaction  of  a  shock  with  an  oblique  interface  between  two  gases.  Re¬ 
sults  of  his  regular-refraction  experiments  agree  closely  with  the  theoretical 
results  of  Polachek  and  Seeger.^  The  experiments  also  revealed  several 
types  of  irregular  refraction  process,  for  which  no  comparison  theory  ex¬ 
isted.  Jahn  discussed  these  patterns  and  showed  that  they  could  be  explained 
qualitatively  by  application  of  simple  gas-dynamic  principles.  His  discussions, 
in  amplified  form,  are  also  given  in  a  series  of  Princeton  University  reports.^ 
The  experimental  set-up  consisted,  ideally,  of  a  two-dimensional  channel 
in  which  there  was,  initially,  an  oblique,  essentially-massless  diaphragm  sep¬ 
arating  two  gases  in  equilibrium.  A  plane,  steady  shock  was  allowed  to  ap¬ 
proach  through  one  of  the  gases  and  interact  with  the  interface,  and  the  re¬ 
sult  was  photographed  during  the  interaction.  As  such,  the  experiments  dif¬ 
fered  from  the  theoretical  model  assumed  by  Polachek  and  Seeger,  and  by 
12 

Taub  in  his  similar  calculations.  In  the  model,  the  interface  was  of  in¬ 
finite  extent,  and  no  effects  from  the  comer  were  considered.  In  the  ex¬ 
perimental  work,  Jahn  was  able  to  separate  the  comer  effects  from  those 
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of  the  pure  refraction  by  suitable  adjustment  of  the  angle  of  the  channel  wall 
just  beyond  its  intersection  with  the  comer.  These  corner  effects  cannot  be 
completely  eliminated,  however,  and  any  complete  theoretical  description  of 
the  processes  must  include  them. 

Because  no  initial  assumptions  concerning  the  nature  of  the  interaction 
are  required  in  PIC-method  calculations,  it  was  expected  that  both  regular 
and  irregular  refraction  processes,  together  with  corner  effects,  could  be 
computed.  Some  representative  calculations  were  carried  out,  and  the  re¬ 
sults  confirmed  the  expectations  and  exhibited  well  certain  properties  of  the 
method. 

In  the  calculations  the  gases  were  considered  to  be  nonviscous,  non¬ 
conducting,  and  polytropic.  They  were  confined  between  rigid  parallel  walls 
in  a  two-dimensional  channel,  with  shock  input  on  the  left  and  a  continuative 
boundary  on  the  right.  Each  was  initially  homogeneous  and  the  two  were  at 
the  same  pressure.  The  section  of  channel  studied  was  divided  into  1200 
square  computational  cells,  having  a  maximum  average  of  four  particles  per 
cell.  The  interface  between  gases  was  inclined  at  45°  in  all  calculations. 
This  was  the  easiest  angle  to  represent  in  the  mesh  of  square  cells;  the 
angle  could  be  changed  if  rectangular  cells  were  allowed,  or  if  the  particle 
placements  were  altered.  With  the  angle  chosen,  however,  all  features  of 
interest  for  this  study  were  revealed. 

The  best  representation  of  an  input  shock  is  obtained  if  the  average 
number  of  particles  per  cell  behind  it  is  an  integer.  Therefore,  in  each 
case  the  density  ratio  across  the  incoming  shock  was  required  to  be  four 
to  two. 

The  results  of  the  calculations  are  presented  mainly  in  the  form  of 
illustrations  of  the  flow  configuration  at  suitable  times  during  the  inter¬ 
actions.  Each  picture  was  formed  by  plotting  the  positions  of  all  the  mass 
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points  in  the  system  and  superimposing  an  Interpretation  of  the  main  features 
by  lines  representing  shocks,  slip  planes,  and  rarefaction  fronts.  These 
positions  were,  in  every  case,  determined  by  observation  of  deflections  of 
the  mass-point  lines.  Where  the  signals  were  weak,  the  results  of  this 
procedure  were  not  always  in  good  agreement  with  known  results;  the  details 
are  discussed  in  the  individual  cases,  and  it  is  pointed  out  that  there  are 
usually  other  satisfactory  means  of  locating  the  weak-signal  positions.  Cer¬ 
tain  peculiarities  of  the  mass-point  plots  are  discussed  below  in  connection 
with  Fig.  E-l. 

Numerous  other  results  were  also  obtained  from  the  computations. 

These  included  internal  and  kinetic  energies  of  each  gas  and  other  functionals 
of  motion  such  as  the  vertical  and  horizontal  momenta.  Because  these  were 
considered  to  be  of  less  interest  here,  they  are  discussed  only  briefly  in  a 
few  cases.  Among  the  least  reliable  results  from  the  calculations  are  such 
quantities  as  local,  instantaneous  densities.  These  and  other  local  features 
which  depend  upon  the  number  of  mass  points  per  cell  may  fluctuate  rather 
strongly  about  their  true  values.  It  is  an  essential  feature  of  the  computing 
method,  however,  that  the  gross  functionals  of  motion  (those  which  extend 
over  numerous  cells  or  depend  upon  time  averaging)  are  well  reproduced  In 
spite  of  these  local  fluctuations. 

In  our  computations,  the  gases  have  been  considered  to  be  nonviscous, 
non-heat-conducting,  and  polytropic.  Effects  from  the  computing-method  ap¬ 
proximations,  however,  can  be  interpreted  as  imparting  certain  other  char¬ 
acteristics  to  the  gases.  Principal  among  these  is  an  effective  "viscosity” 
which  allows  for  the  calculation  of  shocks,  but  which  also  results  in  a  shear 
adhesion.  This  is  mathematically  similar  to  true  viscosity  but  differs,  for 
example,  in  being  anisotropic  in  a  manner  dependent  upon  the  orientation  of 
the  coordinate  system.  The  effects  are  clearly  visible  in  the  results  reported 
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in  this  chapter,  but  are  unimportant  in  their  effect  on  most  features  of 
interest. 

We  have  used  the  notation  of  Jahn,  according  to  which  the  initial  shock 
I  of  strength  £  (£  =  ratio  of  pressure  ahead  to  pressure  behind)  is  incident 
from  the  left  on  the  interface  O.  In  the  resulting  configuration,  there  is  a 
reflected  shock  RS  or  reflected  rarefaction  RR,  a  deflected  interface  D,  and 
a  transmitted  shock  T.  In  Table  II-l  are  tabulated  the  essential  features  of 
the  problems. 


Table  H-l 

LIST  OF  COMPUTATIONS 

Problem  Gas 


Number 

Combination 

PR/PL 

yr 

£ 

1 

Air/  COg 

1.529 

1.405 

1.304 

0.362 

2 

Air/CH 

4 

0.554 

1.405 

1.310 

0.362 

3 

Air/Neon 

0.696 

1.405 

1.667 

0.362 

4 

Air/Polytropic 

Gas 

0.237 

1.405 

1.667 

0.362 

5 

Air/He 

0.138 

1.405 

1.667 

0.362 

B.  Regular  Refraction  with  Corner  Effects 

With  the  restrictions  previously  mentioned,  it  was  not  possible  to  cal¬ 
culate  for  any  situation  corresponding  exactly  to  an  experiment  by  Jahn. 
Instead,  as  a  check  on  our  results  we  calculated  points  on  two  of  the  curves 
of  Polachek  and  Seeger,  choosing  curves  on  which  several  points  had  been 
verified  by  Jahn's  results. 

The  interaction  configuration  for  problem  1  is  shown  in  Fig.  H-l. 
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This  regular  refraction  pattern  consists  of  transmitted  and  reflected  shocks, 
and  a  smoothly  deflected  interface.  Also  visible  is  the  corner  signal  (a 
rarefaction  wave)  which,  in  its  propagation  along  the  interface  through  the 
air,  induced  a  more  slowly  traveling  rarefaction  signal  in  the  CO^.  Inter¬ 
action  of  the  oblique  transmitted  shock  with  the  channel  wall  resulted  in  a 
Mach  reflection  pattern.  The  effects  of  various  weaker  disturbances  can  be 
seen  in  the  pattern  of  mass  points  in  the  lower  part  of  the  air. 

There  are  certain  features  of  the  mass-point  plot  and  interpretation 
which  also  apply  to  the  other  illustrations.  First,  it  is  apparent  that  there 
is  a  vertical  discontinuity  in  the  pattern  of  the  plus  points.  Those  to  the 
right  of  the  discontinuity  were  present  in  the  calculation  region  at  initial 
time,  those  to  the  left  entered  subsequently,  and  the  manner  of  their  intro¬ 
duction  was  not  such  as  to  produce  the  same  pattern  as  resulted  from  com¬ 
pressing  the  gas  already  present.  The  deflections  of  the  vertical  or  diagonal 
mass-point  lines  are  more  significant  in  the  determination  of  boundaries  in 
the  flow  pattern  than  are  those  of  the  horizontal.  Indeed,  near  the  deflected 
interface,  at  which  there  should  be  slippage,  the  effective  "viscosity"  arising 
from  the  computing  method  caused  rather  strong  deflections  of  the  horizontal 
mass-point  lines. 

The  most  clearly  discernible  features  are  the  interface  position  and 
the  positions  of  the  initial  and  transmitted  shocks.  Location  of  the  corner 
and  reflected  signals  is  usually  difficult  and  requires  the  use  of  a  straight¬ 
edge  placed  along  the  mass-point  lines.  In  problem  1,  the  reflected  signal, 
a  shock,  was  extremely  weak.  As  a  result,  its  location  as  determined  by 
the  mass-point-line  deflections  is  not  as  high  as  it  should  be;  the  angular 
discrepancy  is  about  7°.  Reference  to  the  detailed  cellwise  listings  of  the 
calculation,  however,  shows  that  a  signal  had  actually  penetrated  higher  than 
is  shown  in  Fig.  n-1,  and  good  agreement  with  the  correct  result  can  be 
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obtained  by  drawing  a  line  through  the  group  of  uppermost  cells,  each  of 
which  shows  any  change  in  temperature  from  that  in  the  cells  above  it.  The 
positions  of  the  interface  and  of  the  transmitted  shock,  however,  are  in  good 
agreement  with  the  correct  positions,  differing  in  angular  deflection  by  ap¬ 
proximately  2°  and  0°  respectively.  The  position  of  the  comer  signal  is 
likewise  consistent  with  the  expected  position  as  determined  by  the  sound 
speed  in  the  air  behind  the  reflected  shock. 

The  results  for  problem  2  are  shown  in  Fig.  n-2.  In  this  case,  the 
computed  angular  deflection  of  the  interface  was  too  great  by  about  3°,  while 
that  of  the  transmitted  shock  was  too  small  by  about  2°.  Again,  the  reflected 
signal  (a  rarefaction)  was  too  low,  as  measured  by  mass-point-line  deflections, 
with  an  angular  discrepancy  of  about  7°.  The  primary  signals  from  the  cor¬ 
ner  are  more  complicated  in  this  calculation.  They  include  a  weak  signal 
which  propagated  rapidly  through  the  methane  and  a  stronger  signal  whose 
propagation  rate  along  the  interface  through  the  air  was  nearly  the  same  as 
that  observed  in  problem  1.  The  strength  of  this  second  signal  (a  shock) 
was  great  enough  so  that  a  following  rarefaction  was  produced;  this,  in  turn, 
induced  a  rarefaction  in  the  methane.  In  the  lower  part  of  the  region,  there 
are  numerous  signals  whose  effects  on  the  interface  are  shown  by  the  flexures 
thereof.  A  plane  of  considerable  slippage  has  been  drawn  in. 

In  both  problems  1  and  2,  the  accuracy  of  the  calculations  is  consistent 
with  the  resolution  to  be  expected  from  this  size  of  finite-difference  cells. 

The  most  apparent  discrepancies  are  in  the  reflected-signal  positions,  when 
they  are  weak;  but,  in  these  and  other  calculations,  such  weak-signal  loca¬ 
tions  can  be  determined  by  reference  to  detailed  listings  of  cellwise  tem¬ 
peratures. 

C.  Irregular  Refraction 

We  have  investigated,  in  particular,  the  irregular  refraction  pattern 
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arising  when  the  right-hand  gas  has  a  higher  sound  speed  than  that  of  the 
left-hand  gas.  We  were  especially  interested  in  studying  cases  similar  to 
those  in  which  the  experiments  showed  the  deflected  interface  to  be  unstable 
[Ref.  9,  Plate  6,  Fig.  14(c)  and  Plate  12,  Fig.  18]. 

The  irregular  refraction  pattern  and  the  unstable  interface  were  both 
observed  in  the  calculation  of  refraction  from  air  to  helium.  To  study  the 
development  of  these,  we  performed  a  sequence  of  calculations.  The  first 
one— a  regular  refraction  from  air  to  neon— is  shown  in  Fig.  n-3.  The  con¬ 
figuration  is  particularly  simple;  the  corner  signals  were  weak  and  the  de¬ 
flected  interface  is  nearly  straight.  The  reflected  shock  was  extremely 
weak  and  is  not  shown.  In  the  lower  part  of  the  air  there  are  indications 
of  a  slip  plane  and  a  rarefaction.  These  became  progressively  more  evident 
in  the  subsequent  two  calculational  results. 

In  the  next  calculation,  the  density  of  the  right-hand  gas  was  reduced 
to  a  value  intermediate  between  those  of  neon  and  helium,  and  the  tempera¬ 
ture  (hence  also  the  sound  speed)  was  increased  for  initial  equilibrium.  As 
a  result,  the  transmitted-shock  speed  along  the  interface  was  greater  than 
that  of  the  input  shock,  and  detachment  occurred.  The  configuration  is  shown 
in  Fig.  H-4.  The  reflected  shock  is  now  evident,  as  is  its  strong  modifica¬ 
tion  by  the  following  rarefaction.  The  slip  plane  and  rarefaction  in  the  lower 
air  were  much  more  strongly  developed  than  in  problem  3.  The  deflected 
interface  showed  some  instability;  adjacent  to  it,  the  right-hand  gas  was 
turbulent,  while  the  air  remained  fairly  stable. 

The  air-to-helium  refraction  pattern,  presented  in  Fig.  n-5,  shows  a 
strong  development  of  the  irregularity  as  well  as  a  great  instability  of  the 
interface  and  of  the  shocked  helium.  The  irregularly  reflected  shock  inter¬ 
acted  with  the  rarefaction,  and  the  resulting  disturbance  was  somewhat  weak¬ 
ened.  The  lower  slip  plane  and  rarefaction  in  the  air  were  even  more 
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strongly  developed  than  before.  The  result  of  continuing  the  calculation  to 
a  later  time  is  shown  in  Fig.  n-6. 

It  is  not  thought  that  the  structure  of  the  deflected  interface  is  cor¬ 
rectly  represented  in  detail.  Just  as  in  an  experiment,  the  exact  nature  of 
the  structure  depends  on  the  form  of  the  initial  interface  irregularities,  as 
well  as  on  the  nature  of  interactions  with  the  complicated,  small-scale  struc¬ 
ture  of  shocks  and  rarefactions.  In  the  calculation,  the  initial  Irregularities 
of  the  interface  correlate  with  the  nature  of  the  finite-difference  mesh,  and 
many  of  the  small-scale  structures  are  not  resolved  by  that  mesh.  Never¬ 
theless,  it  seems  significant  that  the  computed  refractions  are  stable  in  many 
cases  and  become  unstable  in  circumstances  under  which  the  experimental 
patterns  are  also  unstable. 
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Fig.  II- 1  Interaction  configuration  for  problem 
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Fig.  H-2  Interaction  configuration  for  problem 
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Fig.  II-3  Interaction  configuration  for  problem  3. 
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Fig.  n-4  Interaction  configuration  for  problem 
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Fig.  H-5  Interaction  configuration  for  problem 


Fig.  n-6  Late-time  interaction  configuration  for  problem 


CHAPTER  HI 

SHOCK  PASSAGE  THROUGH  A  DISCONTINUOUSLY  ENLARGED  CHANNEL 

A.  Introduction 

Considerable  attention  has  been  paid  to  the  problem  of  determining 
theoretically  the  changes  in  a  shock  as  it  passes  through  a  channel  of  vari¬ 
able  cross  section.  The  problems  for  weak  shocks  or  gradual  area  changes 
have  been  treated  by  numerous  authors  from  several  points  of  view.  Only 
a  moderate  amount  has  been  written,  however,  about  the  more  difficult  prob¬ 
lems  associated  with  the  passage  of  a  strong  shock  through  a  channel  with 
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rapidly  varying  cross  section.  Whitham  has  devised  an  approximation 

method  which  should  be  useful  in  many  cases.  A  different  approach  has 

14 

been  used  successfully  by  Laporte  for  constricted  channels.  It  does  not 
seem  likely,  nevertheless,  that  any  analytical  treatment  which  might  be  de¬ 
veloped  in  the  near  future  will  be  capable  of  describing  the  entire  flow  field 
in  these  complicated  cases,  but  that  solutions  for  particular  situations  will 
have  to  be  obtained  by  use  of  special  numerical  techniques. 

In  this  chapter,  we  describe  the  application  of  the  PIC  method  to  the 
special  problem  of  determining  the  two-dimensional  flow  of  a  strong  shock 
through  a  discontinuously  enlarged  channel  formed  by  two  rigid  planes  which 
are  parallel  except  at  the  discontinuity.  There,  one  of  the  planes  has  a 
double  right-angle  bend,  doubling  the  channel  width.  This  is  one  element, 
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in  simplified  form,  of  some  of  the  problems  which  arise  in  studies  of  shock 
damage  to  building  interiors  and  in  questions  concerning  establishment  of 
flow  patterns  in  a  shock  tube.  In  this  situation,  development  of  the  flow 
field  with  time  can  be  divided  into  two  phases.  In  the  first,  the  shock,  which 
approached  from  the  narrow  side,  diffracts  about  the  corner.  Behind  it, 
there  is  established  a  complicated  structure  which  includes  an  expansion  fan, 
together  with  several  shocks  and  slip  lines.  The  entire  structure  retains 
geometric  similarity  as  it  expands  linearly  with  time.  In  the  second  phase, 
the  shock  interacts  with  the  wall  of  the  enlarged  channel,  is  further  modified, 
and  the  flow  near  the  comer  eventually  becomes  nearly  steady. 

The  features  of  both  phases  are  discussed  in  some  detail  in  this  chap¬ 
ter.  In  addition,  comparisons  with  experiments  are  presented  for  the  first 
phase  of  development.  For  the  experimental  results,,  we  are  grateful  to 
Dr.  Russell  E.  Duff  who  generously  made  available  unpublished  photographs 
of  the  results  of  experiments  performed  with  nitrogen  at  the  Shock  Tube 
Laboratory,  University  of  Michigan,  in  1950. 

B.  The  Computations 

Two  different  gases  were  studied,  nitrogen  and  helium.  They  were 
considered  to  be  polytropic,  nonviscous,  and  non-heat-conducting,  with  specific 
heat  ratios  y  =  1.404  and  y  =  1.667,  respectively. 

In  each  of  the  computations,  the  Mach  number  behind  the  initial  shock 
was  greater  than  unity  so  that  no  changes  in  flow  occurred  in  the  narrow 
section  of  channel,  in  this  instance  to  the  left.  Thus,  the  region  of  com¬ 
putation  was  chosen  to  cover  a  rectangular  section  of  channel  just  to  the 
right  of  the  discontinuity.  In  the  three  problems  with  nitrogen,  the  region 
was  shorter  along  the  channel  length  in  order  that  the  best  possible  resolu¬ 
tion  could  be  obtained  near  the  comer.  In  the  problem  with  helium,  the 
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region  was  longer  along  the  channel  length  in  order  that  down-channel  effects 
could  be  studied. 

At  the  upper  half  of  the  left  boundary,  the  input  represented  the  steady 
flow  behind  the  initial  shock  which,  at  time  t  =  0,  was  at  the  discontinuity. 

At  the  right,  the  boundary  condition  was  continuative.  The  channel  walls 
were  treated  as  rigid  and  reflective  with  perfect  slippage.  At  initial  time, 
the  gas  in  the  computation  region  was  homogeneous  and  at  rest,  represented 
by  two  mass  points  per  cell.  Units  were  scaled  so  that  a  cell  length  was 
unity,  and  the  undisturbed  gas  was  at  unit  density.  In  the  undisturbed  nitro¬ 
gen,  the  sound  speed  was  unity;  however,  the  shock  in  helium  was  of  infinite 
strength,  and  in  this  case  the  material  speed  behind  it  was  initially  unity. 

The  units  of  time  were  accordingly  chosen. 

C.  Development  of  the  Flow  Patterns  in  Nitrogen 

The  larger  section  of  channel  was  40  units  high,  and  the  region  of 
study  extended  30  units  to  the  right  of  the  discontinuity.  Calculations  were 
performed  for  three  different  strengths  of  the  initial  shock. 

The  configuration  of  mass  points  at  time  t  =  12.593  is  shown  in 
Fig.  m-1  for  the  case  in  which  the  Mach  number  behind  the  shock  was 
M  =  1.008.  The  solid  lines  represent  the  shock  and  rarefaction  fronts  as 
deduced  from  the  deflections  of  the  mass-point  lines.  The  long-dashed 
curves  were  taken  from  a  photograph  by  Duff.  The  short-dashed  curve  is 
theoretical  Prandtl-Meyer  streamline.  The  mass  points  plotted  by  dots  are 
those  which  were  in  the  region  at  t  =  0;  those  plotted  by  pluses  entered 
subsequently.  (The  individual  mass  of  a  plus  point  was  not  the  same  as 
that  of  a  dot  point,  so  that  the  apparent  density  discontinuity  between  them 
is  not  real.) 

A  similar  configuration  is  shown  in  Fig.  m-2  for  time  t  =  6.329.  In 
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this  case,  the  Mach  number  for  flow  behind  the  shock  was  M  =  1.588. 

In  both  cases,  the  shape  of  the  diffracted  shock  is  represented  to 
within  the  expected  resolution  of  the  computations,  while  the  rarefaction 
front,  as  determined  by  deflection  of  the  mass-point  lines,  falls  somewhat 
below  the  experimentally  observed  position.  The  discrepancy  is  not  sur¬ 
prising.  Where  the  upward-traveling  signal  was  weakest,  the  deflections 
of  the  mass  points  were  so  small  as  to  be  undetectable  in  the  plot.  Refer¬ 
ence  to  detailed  cellwise  listings  of  results  from  the  computations  shows, 
however,  that  a  signal  had  indeed  penetrated  higher  than  the  position  shown 
by  the  particle -line  deflections,  and  approximate  agreement  with  the  experi¬ 
mentally  observed  line  can  be  obtained  by  an  appropriate  interpretation  of 
the  listings.  The  situation  is  the  same  as  encountered  with  the  reflected 
signals  in  Chapter  n. 

It  is  not  expected  that  the  Prandtl-Meyer  streamline  should  lie  along 
a  mass-point  line  in  the  early  development  of  the  flow.  The  first  tendency 
of  the  mass-point  lines  is  to  curve  downward;  for  weak  shocks,  this  tendency 
results  in  the  formation  of  the  well-known  spiral  vortex.  As  the  flow  con¬ 
tinues,  however,  the  mass-point  lines  should  approach  the  streamlines.  This 
expectation  was  closely  realized  in  the  two  computations  which  were  extended 
to  late  times.  In  them  (one  case  is  shown  in  Fig.  HI-6),  the  hand-computed 
Prandtl-Meyer  streamlines  lie  very  close  to  the  mass-point  lines  up  to  the 
point  where  the  flow  is  perturbed  by  reflection  from  the  channel  bottom. 

The  configuration  for  M  =  1.296  is  shown  in  a  slightly  different  way 
in  Fig.  m-3.  The  time  is  t  =  10.063.  The  fine  lines  connect  mass-point 
lines  which  formed  a  square  grid  in  the  initially  unshocked  gas,  as  shown 
in  the  lower  right-hand  comer.  The  computed  shock  and  rarefaction  fronts 
were  again  determined  by  the  deflection  positions  of  these  lines.  The  ex¬ 
perimental  positions  are  shown  as  dashed  lines;  also  included  in  this  figure 
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are  the  details  of  additional  flow  structure  observed  on  the  photograph 
(Fig.  HI-4)  near  the  diffracting  comer.  On  the  photograph  the  two  lines 
diverging  from  the  corner  represent  the  end  of  the  expansion  fan  and  the 
slip  stream.  The  former  is  truncated  by  a  shock.  These  structures  are 
suggested  by  the  computations,  but  the  lack  of  resolution  with  this  coarse 
cell  size  precluded  reproduction  of  the  details.  The  structure  just  below 
the  comer  is  not  clearly  correlated  with  any  feature  of  the  computations. 

The  experimentally  observed  curved  line  behind  the  shock,  and  approximately 
parallel  to  it,  is  represented  in  the  computation  by  a  thin  region  of  strong 
slip  and  of  compression  gradient,  as  shown  by  the  mass-point-line  positions. 

Structures  similar  to  these  are  also  present  in  photographs  (not  shown) 
for  the  other  two  Mach  numbers,  corresponding  to  Figs.  HI-1  and  HI-2,  but 
likewise  in  these  cases  the  computation  did  not  well  reproduce  the  details. 

Figures  HI-1  and  m-2  show  a  region  of  instability  in  the  flow  field 
just  below  the  comer.  The  boundary  between  incoming  gas  and  that  which 
was  initially  in  the  computation  region  folds  back  and  forth  in  a  manner 
suggesting  Helmholtz  instability.  Positions  of  this  line  at  two  different  times 
are  shown  in  Fig.  HI-5  for  the  computation  with  M  =  1.296.  This  computa¬ 
tion  was  continued  to  considerably  later  times.  Shortly  after  t  =  20.127,  the 
regularly  folded  structure  of  the  interface  was  no  longer  recognizable,  and 
by  t  =  35  the  region  of  mixing  had  reached  the  right  side  of  the  computation 
region. 

By  t  =  51.178,  considerable  of  the  original  gas  was  still  trapped  in  the 
lower  left  comer.  The  configuration  of  mass  points  at  this  late  time  is 
shown  in  Fig.  HI-6.  The  position  of  the  rarefaction  front  where  it  is  weak 
is  not  well  indicated  by  the  mass-point -line  deflections.  Some  effect  of  the 
reflection  of  the  rarefaction  from  the  upper  channel  wall  is  suggested  by 
the  appearance  of  the  upper  right  comer  of  the  region. 
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D.  Functionals  of  Motion 


Until  a  signal  from  the  corner  has  reached  the  top  or  bottom  of  the 
channel,  all  quantities  in  the  disturbance  region  should  be  functions  of  only 
the  ratios  of  the  position  coordinates  to  the  time,  and  not  otherwise  depend 
on  time.  As  a  consequence,  such  functionals  as  vertical  momentum  and  ex¬ 
cess  kinetic  energy  (over  input)  should  increase  exactly  quadratically  with 
time.  The  first  of  these  can  be  computed  by  hand,  and  the  result  serves 
as  a  check  of  the  machine  computation.  The  result  for  the  case  with 
M  =  1.296  is  shown  in  Fig.  III-7.  In  the  early  stages,  the  computed-momen- 
tum  curve  is  indistinguishable  from  the  theoretical  one,  whose  extension  is 
the  dashed  line.  At  t  =  10.7,  the  shock  arrived  at  the  right-hand  end  of  the 
computation  region  and,  shortly  thereafter,  it  began  to  interact  with  the 
channel  bottom.  Thus,  the  c omputed- momentum  curve  departs  from  the 
theoretical  one.  Much  later,  the  flow  in  the  computation  region  was  in 
nearly  a  steady  state,  at  which  time  the  vertical  momentum  was  nearly  con¬ 
stant. 

In  Fig.  m-7  there  also  is  shown  the  average  specific  internal  energy 
of  a  5  x  5  square  at  the  lower  left  of  the  computation  region.  The  effect 
of  shock  arrival  just  before  t  =  12  is  clearly  visible. 

In  each  of  the  three  problems,  the  total  kinetic  energy  in  the  compu¬ 
tation  region  rose  almost  precisely  at  the  rate  of  the  flux  across  the  input 
boundary  during  the  first  phase  of  development,  indicating  that  kinetic  and 
internal  energy  were  nearly  conserved  separately  during  this  time.  Thus, 
the  increase  in  kinetic  energy  which  the  gas  received  in  the  expansion  fan 
must  have  been  matched  closely  by  the  decrease  as  the  gas  decelerated  at 
the  back  of  the  curved  shock. 
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E.  Infinite  Shock  in  Helium 


A  calculation  similar  to  those  described  above  was  performed  for  an 
infinite  shock  in  helium.  In  this  problem,  the  enlarged  section  of  channel 
was  30  units  high,  and  the  region  of  computation  extended  40  units  along  the 
channel.  The  flow  pattern  for  time  t  =  30  is  shown  in  Fig.  HI-8  and  that 
for  time  t  =  70  in  Fig.  HI-9.  With  the  elongated  calculation  region,  it  was 
possible  to  see  the  reflected  shock  from  the  turbulent  region  along  the  channel 
wall.  Local  mean-velocity  vectors,  whose  magnitudes  can  be  compared  with 
the  input-velocity  vector,  show  a  divided  flow  in  the  turbulent  region.  In  the 
lower  left  comer,  the  fluid  was  slowly  and  irregularly  rotating;  the  upper 
part  interacted  irregularly  with  the  main  stream,  and  occasionally  small 
amounts  were  rapidly  carried  away,  some  being  fed  back  into  the  vortex 
and  the  rest  being  carried  downstream. 
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FLOW  INPUT 


Fig.  in-1  Configuration  of  mass  points  at  time  t  =  12.593  for  the  calcula¬ 
tion  for  nitrogen  with  M  =  1.008.  Solid  and  long-dashed  lines 
represent,  respectively,  the  computed  and  observed  positions  of 
shock  and  rarefaction  fronts.  Short-dashed  line  is  a  theoretical 
Prandtl-Meyer  streamline. 
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FLOW  INPUT 


Fig.  m-2  Configuration  of  mass  points  at  time  t  =  6.329  for  the  calcula¬ 
tion  for  nitrogen  with  M  =  1.588.  Solid  and  long-dashed  lines 
represent,  respectively,  the  computed  and  observed  positions  of 
shock  and  rarefaction  fronts.  Short-dashed  line  is  a  theoretical 
Prandtl-Meyer  streamline. 


Fig.  Ill— 3  Configuration  of  some  of  the  mass-point  lines  at  time  t  =  10.063 
for  the  calculation  for  nitrogen  with  M  =  1.296.  Heavy  solid 
lines  represent  shock  and  rarefaction  front  positions.  Dashed 
lines  represent  these,  together  with  additional  observed  structure, 
as  derived  from  a  photograph,  Fig.  m-4. 


Fig.  m-4  Schlieren  photograph  by  Duff  of  the  flow  field  structure  in  nitro¬ 
gen  with  M  =  1.296. 
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Fig.  m-5  Boundary  between  incoming  gas  and  that  which  was  in  the  region  at  t  =  0  for  the  calcula¬ 
tion  for  nitrogen  with  M  =  1.296.  Only  the  lower  half  of  the  computation  region  is  shown. 


FLOW  INPUT 


(z.oi  x )  wruN3woifti  ay vmnmoq 


Downward  momentum  of  the  nitrogen  in  the  computation  region  for  the  problem  with 
M  =  1.296.  Dashed  spur  is  extension  of  theoretical  curve  for  early  times.  Also  shown 
is  the  average  specific  internal  energy  in  a  5  x  5  square  at  the  lower  left  of  the  com¬ 
putation  region. 


Configuration  of  the  flow  pattern  at  t  =  30  for  the  calculation  for  helium,  showing  mass 
point  lines,  shock  front,  Mach  line  from  corner,  and  boundary  of  turbulent  region. 


Configuration  of  the  flow  pattern  at  t  =  70  for  the  calculation  for  helium,  showing  mass- 
point  lines,  region  of  turbulence  with  shock  reflected  therefrom,  Mach  line  from  corner, 
and  some  local -velocity  vectors. 


CHAPTER  IV 

INTERACTION  OF  A  SHOCK  WITH  A  DEFORMABLE  OBJECT 

A.  Introduction 

The  encounter  of  a  strong  shock  with  a  deformable  object  can  result 
in  a  very  complicated  flow  field  which  contains  both  diffraction  and  refraction 
processes.  The  PIC  method  calculations  have  been  applied  to  a  study  of 
several  idealized  situations  in  which  the  object  and  the  surrounding  gas  were 
both  represented  by  the  polytropic  equation  of  state  of  a  monatomic,  non- 
viscous,  nonconducting  gas.  The  flow  was  two-dimensional,  confined  between 
rigid  parallel  channel  walls.  The  object  was  rectangular  in  shape  and  at¬ 
tached  to  one  of  the  walls;  It  and  the  surrounding  gas  were  initially  cold, 
so  that  the  incoming  shock  strength  was  infinite. 

In  all  calculations,  units  were  scaled  in  such  a  way  that  the  width  of 
a  computational  cell  was  unit  distance,  and  the  material  velocity  behind  the 
input  shock  was  unity.  The  initial  density  of  gas  around  the  object  was  like¬ 
wise  unit  mass  per  cell,  while  the  object  density  was  four  times  as  great. 
Thus,  the  speed  of  the  incoming  shock  was  1.33  and  the  density  behind  it 
was  4.  The  channel  was  24  cells  wide  and  the  computation  region  in  it  was 
50  cells  long.  The  object  was  10  cells  wide  in  each  case,  but  its  length  was 
variable.  Along  the  left  boundary  there  was  input  corresponding  to  constant 
conditions  behind  the  shock. 
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B.  Configurations  of  the  Flow  Fields 

A  typical  interaction  configuration  is  shown  in  Fig.  IV-1.  The  object 
was  a  long  rectangle,  extending  to  the  right  side  of  the  computation  region. 
Positions  of  the  interface  and  of  the  shocks  were  determined  by  reference 
to  plots  of  the  mass  points,  in  a  manner  similar  to  that  used  in  previous 
chapters. 

A  similar  result  is  shown  in  Fig.  IV-2.  In  this  case  the  object  length 
was  twice  its  width.  In  addition  to  the  features  shown  in  Fig.  IV-1,  there 
is  also  shown  the  set  of  initially-horizontal  mass-point  lines.  These  are 
dashed  in  the  turbulent-vortex  region,  where  the  flow  pattern,  as  represented 
by  the  particle  positions,  is  considerably  more  contorted. 

C.  Functionals  of  Motion 

In  the  various  computations  reported  in  this  paper,  numerous  functionals 
of  the  motion  were  calculated.  For  this  set  of  calculations,  it  is  appropriate 
that  they  be  discussed  in  some  detail,  because  they  demonstrate  several  prop¬ 
erties  of  the  PIC  method.  In  particular,  several  of  the  functionals  can  be 
compared  with  corresponding  hand-computation  results.  Thus,  for  the  prob¬ 
lem  whose  late-time  configuration  is  shown  in  Fig.  IV-1,  the  vertical  momen¬ 
tum  of  the  entire  flow  field  could  be  computed  exactly  by  hand  for  times  up 
to  when  a  signal  from  the  comer  of  the  object  reached  one  of  the  channel 
walls.  The  comparison  is  shown  in  Fig.  IV-3.  The  effect  of  collision  of  a 
comer  signal  with  the  channel  wall  is  clearly  visible  at  t  =  26. 

Likewise,  the  circulation  around  the  computation  region  could  be  hand- 
computed  as  a  function  of  time.  This  was  done  for  the  shorter  object,  since 
the  result  shows  a  discontinuity  in  slope  resulting  from  shock  break-through 
at  the  back  edge  just  after  t  =  28.  The  result  is  also  shown  in  Fig.  IV-3. 

The  change  in  slope  is  clearly  visible  in  the  machine-computed  results. 
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Finally,  an  approximate  hand  computation  could  be  made  of  the  kinetic 
and  internal  energies  of  the  object;  the  results  should  be  valid  for  early 
times  after  the  encounter.  The  comparison  is  shown  in  Fig.  IV-4  for  the 
short-object  calculation.  The  kinetic  energy  from  the  machine  computation 
agrees  with  the  hand-computed  result  nearly  as  well  as  do  the  circulation 
and  vertical  momentum.  The  time-wise  lag  in  the  internal  energy  curve 
has  also  been  observed  in  various  other  PIC  method  calculations.  It  is 
easily  explained  by  the  nature  of  the  finite-difference  equations,  and  can  be 
made  smaller  by  a  decrease  In  cell  size. 
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Fig.  IV-1  Interface  and  shock  positions  at  time  t  =  35  for  the  long  object. 


CIRCULATION 


Fig.  IV-3  Vertical  momentum  and  circulation  of  the  entire  computation 

region  as  functions  of  time  for  the  long-  and  short-object  calcu¬ 
lations,  respectively.  Datum  points  are  from  the  machine  com¬ 
putation;  solid  curves  were  hand-computed,  using  shock  theory. 
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VERTICAL  MOMENTUM 


Fig.  IV-4  Internal  and  kinetic  energies  of  the  object  as  functions  of  time 
for  the  short-object  calculation.  Datum  points  are  from  the 
machine  computation;  solid  curve  was  hand-computed  and  repre¬ 
sents  both  internal  and  kinetic  energies. 
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CHAPTER  V 

HYPERSONIC  SHEAR  FLOW  WITH  PERTURBED  INTERFACE 


A.  Introduction 

We  have  used  the  PIC  method  to  study  the  plane,  two-dimensional  in¬ 
teraction  of  two  gases  moving  past  a  perturbed  slip  plane.  The  gases  were 
confined  between  infinite,  parallel,  rigid  walls;  in  cross  section,  the  initial 
slip  plane  was  approximately  a  low-amplitude  sine  wave  with  mean  position 
halfway  between  the  walls.  The  gases  were  initially  cold  (zero  sound  speed); 
the  upper  one  was  moving  to  the  right  and  the  lower  one  to  the  left,  both 
parallel  to  the  channel  walls  and  with  the  same  speed.  Both  gases  were 
polytropic,  nonviscous,  nonconducting,  and  monatomic  (specific  heat  ratio 
y  =  5/3). 

The  initial  configuration  was  perfectly  periodic  along  the  channel,  and 
it  was  assumed  that  subsequent  interactions  would  preserve  that  periodicity. 
One  period  was  divided  into  1200  square  cells— 30  from  wall  to  wall  and  40 
along  the  channel.  The  total  number  of  mass  points  was  4800.  Boundary 
conditions  were  reflective  at  the  rigid  channel  walls  and  periodic  at  the  ends. 

Dimensions  were  scaled  so  that  each  cell  was  of  unit  width,  the  initial 
velocities  were  of  unit  magnitude,  and  the  density  of  the  upper  gas  was  four 
units.  Thus,  the  time  unit  was  the  time  required  for  the  undisturbed  gas  to 
move  one  cell.  The  initial  perturbation  of  the  slip  line  was  always  the 
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same— of  unit  amplitude  and  in  the  form  of  a  step-function  approximation 
(along  cell  boundaries)  to  a  sine  wave  with  origin  at  the  left  of  the  period. 
The  parameter  which  was  varied  among  calculations  was  R,  the  ratio  of 
density  of  the  upper  to  lower  gases.  The  values  studied  were  R  =  1,  0.5, 
and  0.2,  and  the  scaling  laws  allowed  additional  results  to  be  derived  from 
these,  equivalent  to  R  =  2  and  5. 

B.  The  Interaction 

Along  the  section  of  slip  line  with  positive  slope,  the  gases  collided 
in  the  early  stages.  From  this  collision  line,  shocks  proceeded  into  each  gas; 
the  gases  were  heated,  and  vertical  kinetic  energy  was  created.  The  shocks 
diverged  and  were  carried  along  as  shock  pulses  in  undulating  ribbon  form, 
separated  by  rarefied  regions  of  relative  stagnation.  In  the  early  stages,  a 
cavity  opened  along  the  section  of  slip  line  with  negative  slope. 

Typical  appearance  of  a  well-developed  configuration  is  shown  in 
Fig.  V-l  for  the  initial  density  ratio  R  =  1;  the  elapsed  time  was  30.7  units 
(the  free-stream  motion  had  carried  the  gas  about  three-fourths  of  a  period). 
Mass-point  lines  shown  for  the  upper  gas  were  initially  horizontal  and  evenly 
spaced;  those  for  the  lower  gas  were  initially  vertical  (transverse  to  flow) 
and  also  evenly  spaced.  The  dashed  line  shows  the  mean  Interface  between 
the  gases;  initially,  this  approximately  followed  the  dotted  line.  The  central 
strip  is  rarefied  and  strongly  turbulent  as  indicated  by  the  admedian  mass- 
point  line  of  the  upper  gas.  For  R  =  1,  the  configuration  always  satisfied 
the  symmetry  property  that  vectors  reflected  through  the  center  point  should 
change  sign  but  not  magnitude,  while  scalars  remain  unchanged. 

In  Fig.  V-2  are  shown  the  positions  of  the  diverging  shock  fronts  at 
various  times  for  the  problem  with  density  ratio  R  =0.2.  These  locations 
were  defined  by  the  positions  at  which  the  magnitude  of  vertical  velocity 
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was  one-tenth  of  the  Initial  free-stream  speed.  The  maximum  and  minimum 
heights  of  the  two  shock  fronts  are  shown  as  functions  of  time  in  Fig.  V-3. 
By  t  =  30,  these  shock  fronts  were  moving  at  speeds  of  0.36  and  0.14  cells 
per  time  unit  in  the  upper  and  lower  gases,  respectively. 

In  addition  to  the  two  diverging  shock-pulse  ribbons,  there  was  a  cen¬ 
tral  region  of  rarefaction  wherein  both  gases  had  dropped  to  about  one-third 
normal  density  at  late  times.  In  Fig.  V-4  is  shown  the  compression  as  a 
function  of  height  above  lower  channel  wall,  at  time  t  =  30,  for  the  problem 
with  density  ratio  R  =  0.2.  These  compressions  were  obtained  in  the  com¬ 
putation  as  cellwise  quantities,  "quantized"  by  integral  numbers  of  particles 
in  each  cell,  and  have  been  smoothed  in  the  plots. 

The  distinction  between  the  diverging  ribbons  and  the  central  rarefaction 
region  is  also  strongly  indicated  by  the  vertical  profiles  of  specific  internal 
energy.  One  of  these  is  shown  in  Fig.  V-5  at  time  t  =  30,  at  channel  mid¬ 
length,  for  the  density  ratio  R  =  0.2  (corresponding  to  the  upper  part  of 
Fig.  V-4).  The  ribbon  temperature  was  slightly  higher  than  that  behind  a 
theoretical  plane  shock  of  this  speed.  (A  shock  speed  of  0.36  would  indicate, 
with  y  =  5/3,  a  vertical  material  speed  of  0.27— which  is  close  to  the  ob¬ 
served  value— but  this  in  turn  leads  to  a  specific  internal  energy  of  0.37.) 

The  explanation  lies  in  the  fact  that  the  shock  was  actually  oblique  on  its 
left  face,  where  the  shock  speed  was  greater  relative  to  the  material. 

Energy  was  transferred  from  the  denser  to  lighter  gas  at  a  rate  which 
is  roughly  proportional  to  (1  —  R)/R  in  the  range  of  density  ratios  considered. 
This  transferal  was  a  secondary  effect,  in  that  the  rate  remained  negligible 
for  some  time  and  then  gradually  increased. 

Another  secondary  energy  effect  was  the  production  of  vertical  kinetic 
energy.  The  rate  of  production  roughly  equaled  the  rate  at  which  energy 
was  transferred  to  the  lighter  gas  in  the  case  R  =  0.2,  but  was  somewhat 
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larger  for  larger  values  of  R.  Crude  calculations  using  simple  shock  theory 
suggest  that  the  rate  of  production  of  vertical  kinetic  energy  in  the  upper 

_3 

gas  should  be  proportional  to  (1  +  VR)  for  all  values  of  R.  This  crude 
calculation  agrees  surprisingly  well  with  the  results  of  the  full  machine  cal¬ 
culations  for  early  times.  In  Fig.  V-6  is  shown  the  vertical  kinetic  energy 
in  the  upper  gas  for  five  different  densities  of  the  lower  gas.  (The  curves 
for  R  =  2.0  and  5.0  were  obtained  by  scaling  the  energies  of  the  lower  gas 

from  the  runs  with  R  =  0.5  and  0.2,  respectively.)  Also  shown  as  a  set  of 

-3 

isolated  points  at  t  =  10  are  the  values  of  45(1  +  VR)  ,  the  constant  of 
proportionality  having  been  chosen  to  fit  the  value  for  R  =  0.2.  Variation 
in  the  accuracy  of  agreement  at  earlier  times  is  consistent  with  the  expected 
variations  in  the  finite -difference  results.  The  late-time  drop  in  the  curve 
for  R  =  0.5  arises  from  collision  of  the  shock  with  the  upper  channel  wall. 
The  other  problems  were  not  run  far  enough  for  collision  because  of  the 
machine  time  involved. 

A  primary  energy  effect  was  the  production  of  internal  energy.  In 
Fig.  V-7  is  shown  the  internal  energy  of  the  upper  gas  as  a  function  of 
time.  The  inset,  at  the  same  scale,  shows  the  early-time  section  of  the 
curves  as  the  data  came  from  the  machine.  The  contortions  can  be  traced 
to  a  fictitious  effect  of  the  finite -difference  technique;  these  have  been 
smoothed  slightly  in  the  full  curves. 
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mass-point  lines;  lower  gas  shows  transverse  mass-point  lines.  Dashed  line  is  interface 
position  which  was  originally  along  the  dotted  line.  Locus  of  ends  of  lower  mass-point 
lines  delineates  normal  density  isopycnic. 


INITIAL  INTERFACE 


/ 


oo  o 

—  CM  ro 

■  mm 


-  73  - 


Position  of  the  maximum  and  minimum  shock  heights  as  a  function  of  time  for  the  prob 
lem  with  R  =  0.2. 
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Fig.  V-5  Vertical  internal  energy  profile  at  channel  midlength  at  time  30  units  for  the  problem 
with  R  =  0.2. 
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Fig.  V-6 


Internal  energy  of  the  upper  gas  as  a  function  of  time  for  various  values  of  the  density 
ratio.  Inset,  at  the  same  scale,  shows  the  unsmoothed  appearance  at  early  times. 


CHAPTER  VI 

TAYLOR  INSTABILITY 

The  PIC  method  has  been  applied  to  problems  involving  the  instability 
of  irregular  gaseous  interfaces  subjected  to  normal  acceleration.  A  typical 
sequence  of  configurations  is  shown  in  Fig.  VI-1.  The  calculations  were  per¬ 
formed  in  cylindrical  coordinates,  with  the  axis  forming  the  left  boundary 
of  the  picture;  the  other  boundaries  were  rigid  and  reflective.  The  two 
gases  were  initially  cold,  nonviscous,  nonconducting  and  polytropic,  with 
specific  heat  ratio  y  =  2.  The  initial  boundary  between  them,  as  shown  in 
Fig.  VI-1,  was  perturbed  to  a  square-tooth  shape.  Smaller  perturbations 
were  not  studied  because  of  lack  of  resolution;  larger  ones  could  not  be 
studied  because  of  the  effects  of  boundary  signals.  Acceleration  was  sup¬ 
plied  by  an  applied  pressure  in  the  empty  cells;  the  procedure  is  discussed 
in  Chapter  I,  Section  B-6.  As  mentioned  there,  the  "empty"  region  behaves 
as  though  it  were  filled  by  a  very  hot  gas  at  very  low  density.  The  pressure 
is  finite  (specified)  but  the  sound  speed  is  very  high,  so  that  the  pressure 
remains  homogeneous.  R  is  assumed  that  this  "gas"  is  backed  by  an  in¬ 
finite  reservoir,  so  that  its  pressure  is  constant  in  time.  Thus  it  is  ex¬ 
pected  that  the  applied-pressure  boundary  will  always  be  unstable. 

A  shock  passing  through  the  upper  gas  eventually  crosses  the  boundary 
between  gases.  The  resulting  compression  in  the  vicinity  of  the  boundary 
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(by  a  factor  of  three  for  y  =  2)  initially  decreases  the  perturbation  amplitude. 
If  the  ratio,  R,  of  density  of  lower  gas  to  that  of  upper  gas  is  greater  than 
unity,  then  it  is  expected  that  the  perturbation  amplitude  will  subsequently 
grow.  For  R  <  1,  the  surface  should  be  stable. 

In  each  of  the  calculations,  there  were  1200  square  computational  cells; 
units  were  scaled  so  that  each  was  of  unit  width.  The  applied  pressure  was 
such  as  to  give  a  material  velocity  behind  the  initial  shock  of  0.1  cell  per 
unit  time. 

In  Fig.  VI-1  there  is  shown  a  sequence  of  configurations  of  the  upper 
gas  for  the  calculation  with  density  ratio  R  =  2.  Corresponding  configura¬ 
tions  are  shown  in  Fig.  VI-2  for  density  ratio  R  =  0.5.  In  the  latter  case, 
the  lower  surface  was  stable,  but  had  changed  phase  by  late  time.  In  the 
former  calculation  the  lower-surface  perturbation  amplitude  had  increased 
back  to  its  original  amount  by  time  t  =  230.  Subsequently,  however,  the 
amplitude  remained  nearly  constant,  while  the  shape  of  the  interface  changed. 
In  both  cases,  the  upper  surface  had  become  perturbed;  the  amplitude  in¬ 
creased  most  in  the  calculation  with  R  =  0.5,  in  which  there  was  the  greatest 
acceleration  of  that  surface. 

Results  of  a  similar,  but  more  extreme,  pair  of  calculations  are  shown 
in  Figs.  VI-3  and  VI-4.  The  initial  configuration  was  the  same  as  before 
in  both  cases,  but  is  not  shown  in  Fig.  VI-3  because  of  overlap  in  drawing 
the  later  configuration.  In  the  first  of  this  pair,  the  density  ratio  was  R  =  20; 
by  time  t  =  230,  the  lower-surface  amplitude  had  increased  to  three-halves 
its  original  amount,  and  the  upper  surface  was  considerably  irregular.  For 
the  case  R  =  0.05,  there  was  again  a  phase  change  in  the  lower-surface 
perturbation.  The  shape  of  the  surface  was  quite  distinctive.  The  upper 
gas  remained  quite  thick;  the  upper  surface  is  not  shown,  again  because  of 
overlap  with  earlier  lower-surface  positions,  and  because  only  by  the  last 
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time  had  that  surface  become  appreciably  perturbed. 

In  one  additional  calculation,  the  lower  gas  was  omitted  and  the  per¬ 
turbation  placed  along  the  upper  surface  of  the  upper  gas,  adjacent  to  the 
applied  pressure.  A  sequence  of  interface  configurations  is  shown  in  Fig. 
VI-5.  Some  small  droplets  which  broke  from  the  top  of  the  axial  spike  are 
not  shown.  The  lowest  configuration  in  the  drawing  is  displaced  downward 
from  its  true  position  to  prevent  overlap. 

It  is  unfortunate  that  none  of  these  results  could  be  compared  with 
experiment  or  other  theory.  Qualitatively,  the  results  appear  reasonable; 
it  is  thought  that  they  also  have  quantitative  accuracy,  to  some  extent. 
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Fig.  VI- 1  Sequence  of  configurations  of  the  upper  gas  for  calculation  with 
density  ratio  R  =  2. 
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Fig.  VI-2  Two  configurations  of  the  upper  gas  for  the  calculation  with 
density  ratio  R  =  0.5. 
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Fig.  VI-5  Sequence  of  configurations  of  the  upper  gas  surface  for  the  cal¬ 
culation  with  applied-pressure-interface  perturbation.  Lowest 
configuration  is  displaced  as  shown  to  avoid  overlap. 
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CHAPTER  VH 


VISCOUS-FLOW  CALCULATIONS 


A.  Introduction 

The  procedure  outlined  in  Chapter  I,  Sec.  C,  has  been  applied  to  sev¬ 
eral  problems  concerning  the  flow  of  a  viscous,  heat-conducting,  polytropic 
gas.  Results  of  several  of  the  calculations  are  presented  here. 

The  gas  in  each  case  was  air  for  which  we  took  y  =  1.4,  A  =  —2/3, 

B  =  1.9.  The  first  viscosity  coefficient  was  assumed  to  vary  with  specific 
internal  energy  (hence  with  temperature)  according  to  the  relation  fx  = 
where  the  exponent  was  a  constant.  (In  Sec.  B  below,  n  =  0.5;  in  Sec.  C, 
n  -  0.)  In  each  case,  the  channel  wall  was  allowed  to  be  conducting  and 
held  at  a  fixed  temperature.  There  were  ten  cells  across  the  channel  and 
three  columns  of  cells  along  one  period  of  the  channel.  In  each  case,  the 
flow  was  actually  one-dimensional,  but  results  were  improved  by  along- 
channel  averaging. 

B,  Couette  Flow 

Initially,  the  gas  was  at  rest  and  at  uniform  density,  p^,  and  specific 
internal  energy,  IQ  =  1.0.  (The  value  of  I  at  the  walls  was  held  fixed  at 
that  value.)  At  time  t  -  0,  the  lower  wall  (at  y  =  0)  was  impulsively  ac¬ 
celerated  to  velocity  u^,  and  thereafter  moved  at  that  constant  rate. 
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It  can  be  shown  that  under  these  circumstances  the  final  (steady-state) 
maximum  value  of  I  should  be 

2 

I  =  1  +^_ 

m  0  8B 

Thus,  we  may  distinguish  between  low-velocity  and  high-velocity  Couette  flow 
according  as  uQ  is  distinctly  less  than  or  greater  than  8BI^.  In  the  former 
case,  the  temperature  (hence  the  density)  remains  nearly  uniform  across  the 
channel.  In  the  latter,  the  central  temperature  is  large  and  most  of  the 
mass  is  confined  to  narrow  bands  at  the  walls. 

At  first  we  placed  the  particles  in  an  orderly  fashion,  four  per  cell, 
one  at  the  center  of  each  quadrant.  In  the  low- velocity  calculations,  the 
results  were  satisfactory  because  the  density  changes  were  so  small  that 
there  were  no  vertical  particle  crossings.  In  the  high-velocity  calculations, 
however,  such  functionals  of  motion  as  total  kinetic  energy  of  the  system, 
plotted  as  a  function  of  time,  showed  discontinuities  in  slope.  These  oc¬ 
curred  whenever  a  row  of  particles  all  simultaneously  made  a  vertical 
crossing.  An  alternate  procedure  in  particle  placement  was  found  to  be 
considerably  more  satisfactory:  Initially,  the  positions  of  the  four  particles 
in  each  cell  were  generated  at  random.  As  a  result,  the  functionals  of 
motion  were  much  smoother  and  in  better  agreement  with  analytical  calcu¬ 
lations  . 

Simple  tests  of  the  calculation  procedure  produced  results  whose  ac¬ 
curacy  is  as  good  as  expected  with  the  coarse  mesh  used.  Some  plots  of 
velocity  and  specific  internal  energy  are  shown  as  functions  of  height  above 
lower  channel  wall  in  Figs.  VH-1  and  VII-2.  This  was  a  low-velocity  case; 

u  =1.0.  The  results  of  a  high-velocity  Couette  flow  calculation  (u  =  10) 

0  u 
are  shown  in  Fig.  VII-3  for  a  late  time.  The  five-times-compression  curve 


-  88 


should  rise  to  a  value  of  36  at  each  side  to  maintain  constant  pressure.  It 
was  apparent  that  the  channel  center  would  never  cool  to  the  theoretical 
limit  in  this  calculation.  The  reason  could  be  traced  to  a  fault  of  the  bound¬ 
ary  conditions  at  the  walls.  It  was  assumed  that  there  the  specific-internal- 
energy  gradient  was  the  value  obtained  by  subtracting  from  the  value  in  the 
first  cell  the  assigned  wall  value  and  dividing  the  result  by  half  the  cell 
width.  The  gradient  thus  calculated  was  much  less  than  the  true  value; 
thus,  a  quadratic  fit  procedure  was  tried  for  calculating  the  gradients.  The 
result  showed  considerably  better  agreement  in  the  value  of  maximum  tem¬ 
perature. 

Even  though  the  final  density  was  far  from  uniform,  the  along-channel 
momentum  from  the  machine  calculation  asymptotically  approached  a  value 
which  was  very  close  to  that  obtained  analytically  as  a  result  of  the  simple 
assumption  that  the  density  and  velocity  gradients  were  constant  across  the 
channel. 

C.  Poiseuille  Flow 

The  initial  conditions  were  the  same  as  in  Sec.  B.  At  time  t  =  0,  a 
body  acceleration  was  applied:  g  =  0.04.  An  approximate  analytical  solution 
was  obtained  for  the  problem,  and  comparisons  with  the  machine -computed 
results  are  shown  in  Figs.  VII-4,  vn-5,  and  VII-6.  The  discrepancies  are 
of  the  same  magnitude  as  the  expected  error  in  the  approximate  solution. 
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HEIGHT  ABOVE  LOWER  WALL 


Fig.  vn-1  Solid  curves  show  machine-calculated  velocity  as  a  function  of 
height  for  low-velocity  Couette  flow.  Dashed  curves  are  ana¬ 
lytical  solutions. 
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SPECIFIC  INTERNAL  ENERGY 


HEIGHT  ABOVE  LOWER  WALL 


Fig.  vn-2  Solid  curves  show  machine-calculated  specific  internal  energy 
as  a  function  of  height  for  low-velocity  Couette  flow.  Dashed 
curve  is  analytical  steady-state  solution. 
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Fig,  VII-3  Specific  internal  energy,  velocity,  and  five  times  the  compres¬ 
sion  plotted  as  functions  of  height  for  a  high-velocity  Couette 
flow  calculation.  Dashed  curve  is  analytical  steady-state  solu¬ 
tion.  Points  show  the  actual  cellwise  values  of  five  times  the 
compression. 
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HEIGHT  ABOVE  LOWER  WALL 

I  Velocity  as  a  function  of  height  in  the  Poiseuil 
Solid  curve  is  from  machine  calculation;  dashei 
annroximate  solution. 


HEIGHT  ABOVE  LOWER  WALL 


Fig.  vn-5  Specific  internal  energy  as  a  function  of  height  in  the  Poiseuille 
flow  problem.  Solid  curve  is  from  machine  calculation;  dashed 
is  analytical  approximate  solution. 


rH-6  Horizontal  momentum  as  a  function  of  time  in  the  Poiseuille-flow  problem.  Solid  curve 
is  from  machine  calculation;  datum  points  are  analytical  approximate  solution. 
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